function [ U, s, V, output ] = IRNNMatCompTest( O, lambda, theta, para )

if(isfield(para, 'maxTime'))
    maxTime = para.maxTime;
else
    maxTime = 1e+9;
end

if(isfield(para, 'decay'))
    decay = para.decay;
else
    decay = 0.96;
end

Omega = sparse(O ~= 0);

maxIter = para.maxIter;
tol = para.tol;
tau = para.tau;
regType = para.regType;

X = full(O);
s = svd(X);
lambdaMax = max(s(:));
s = 1e+6*ones(length(s), 1);
maxR = 5;
R = randn(size(O, 2), maxR);
s = s(1:size(R, 2));

flagTime = tic;
obj = zeros(maxIter, 1);
for i = 1:maxIter
    switch(para.regType)
        case 1
            lambdai = abs(lambdaMax - lambda)*(decay^i) + lambda;
            thetai = theta + (decay^i)*lambdai;
        case 2
            lambdai = lambda;
            thetai = theta;
        case 3
            lambdai = abs(lambdaMax - lambda)*(decay^i) + lambda;
            thetai = theta;
    end
    
    % get weight
    [ w ] = getWeight( s, thetai, lambdai, regType );
    
    Z = X - (1/tau)*(X.*Omega - O.*Omega);
    
    R = randn(size(O, 2), maxR);
    Q = powerMethod(Z, R, 1, 1e-6);
    [U, s, V] = svd(Q'*Z, 'econ');
    s = diag(s);
    s = max(s - w/tau, 0);
    nnzS = sum(s > 1e-6);
    
    % R = V;
    
    U = Q*U(:, 1:nnzS);
    V = V(:, 1:nnzS);
    X = U*diag(s(1:nnzS))*V';
    
    obj(i) = (1/2)*sumsqr((X.*Omega - O.*Omega));
    obj(i) = obj(i) + funRegC(s, length(s), lambda, theta, regType);
    
    if(i > 1)
        delta = obj(i) - obj(i-1);
    else
        delta = inf;
    end
    
    fprintf('iter %d, (obj:%.3d, tol:%.3d), rank %d, lambda %.2d \n', ...
        i, obj(i), -delta, nnz(s), lambdai);
    runTime = toc(flagTime);
    if(abs(delta) < tol || runTime > maxTime)
        break;
    end
end
obj = obj(1:i);

s = s(s > 0);
s = diag(s);
output.obj = obj;
output.rank = nnz(s);

end



